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Abstract Ensemble modeling of coronal mass ejections (CMEs) provides a probabilistic forecast 
of CME arrival time which includes an estimation of arrival time uncertainty from the spread and 
distribution of predictions and forecast confidence in the likelihood of CME arrival. The real-time 
ensemble modeling of CME propagation uses the Wang-Sheeley-Arge (WSA)-ENLIL-I-Cone model 
installed at the Community Coordinated Modeling Center (CCMC) and executed in real-time at 
the CCMC/Space Weather Research Center. The current implementation of this ensemble modeling 
method evaluates the sensitivity of WSA-ENLIL-bCone model simulations of CME propagation to 
initial CME parameters. We discuss the results of real-time ensemble simulations for a total of 
35 CME events which occurred between January 2013 - July 2014. For the 17 events where the 
CME was predicted to arrive at Earth, the mean absolute arrival time prediction error was 12.3 
hours, which is comparable to the errors reported in other studies. For predictions of CME arrival 
at Earth the correct rejection rate is 62%, the false-alarm rate is 38%, the correct alarm ratio is 
77%, and false alarm ratio is 23%. The arrival time was within the range of the ensemble arrival 
predictions for 8 out of 17 events. The Brier Score for CME arrival predictions is 0.15 (where a 
score of 0 on a range of 0 to 1 is a perfect forecast), which indicates that on average, the predicted 
probability, or likelihood, of CME arrival is fairly accurate. The reliability of ensemble CME arrival 
predictions is heavily dependent on the initial distribution of CME input parameters (e.g. speed, 
direction, and width), particularly the median and spread. Preliminary analysis of the probabilistic 
forecasts suggests undervariability, indicating that these ensembles do not sample a wide enough 
spread in CME input parameters. Prediction errors can also arise from ambient model parameters, 
the accuracy of the solar wind background derived from coronal maps, or other model limitations. 
Finally, predictions of the Kp geomagnetic index differ from observed values by less than one for 11 
out of 17 of the ensembles and Kp prediction errors computed from the mean predicted Kp show 
a mean absolute error of 1.3. 
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1. Introduction 


Ensemble modeling has been employed in weather forecasting in order to quantify prediction uncer¬ 
tainties and determine forecast confidence dSivillo, Ahlquist, and Toth, 1997| . Individual forecasts 
which constitute an ensemble forecast represent possible scenarios which approximate a probabil¬ 
ity distribution that reflects forecasting uncertainties. Such uncertainties which when considered 
as a group include those associated with initial conditions (such as observational uncertainties), 
techniques and models. Different forecasts in the ensemble can start from different initial condi¬ 
tions and/or be based on different forecasting models/procedures. In the simplest application, the 
ensemble mean or a weighted mean can be taken as a single forecast. The ensemble mean should 
perform better than individual ensemble members by emphasizing systematic features found in 
all members. However, an ensemble also contains additional information about possible scenarios 
and their probabilities and thus provides a probabilistic forecast. For example, ensemble modeling 
provides a quantitative description of the forecast probability that an event will occur by giving 
event occurrence predictions as a percentage of ensemble size. This conveys the level of uncertainty 
in a given forecast in contrast to a categorical yes/no forecast. Additionally, all ensemble forecast 
members can be plotted together to allow visualization of the uncertainty among ensemble members, 
and their clustering distribution. An example of such a visualization is hurricane track “plume” 
maps in weather forecasting. Regions where members tend to coincide/cluster can be taken to have 
a higher forecast confidence. 

To understand the uncertainties in space weather forecasting, ensemble coronal mass ejection 
(CME) forecasting efforts have now begun in space weather models of the heliosphere. [Fry et aL\ 

(I2003|) , IMcKcnna-Lawlor et al. \ (|2006p , and ISmith et al. \ (|2009p compared the performance of real¬ 
time shock arrival time forecasts following solar events (since 1997) from the three “Fearless 
Forecast” models: Shock Time of Arrival (STOA)(|Dryer, 1974^, Interplanetary Shock Propagation 
Model (ISPM) ( [Smith and Dryer, 1990[ ), and Hakamada-Akasofu-Fry (HAFv.2)( |Dryer et al, 20M ). 

While there are many models predicting the evolution of CMEs (see ( |Zhao and Dryer, 2014[ ) and ref¬ 
erences therein), only the Wang-Sheeley-Arge (WSA) coronal model ( |Arge and Pizzo, 2000||Arge et al., 2004[ )| 
coupled with the global heliospheric ENLIL solar wind model dOdstrcil, 2003 1 has been used exten¬ 
sively in space weather operations world-wide. The first effort in utilizing this model for ensemble 
forecasting of CME propagation was reported bv IPnlkkinen et aPpOlljl . lEmmons et i7i.ip013[l per¬ 
formed WSA-ENLIL ensemble CME modeling using 100 ensemble members for 15 historical events 
with automatically determined cone model CME parameters (Pulkkinen, Oates, and Taktakishvili, 2010().[ 


They found that the observed CME arrival was within the ensemble prediction spread for 8 out 
of the 15 events. ILee et aP (|2013|) discuss ensemble modeling of CME propagation with WSA- 
ENLIL for an event study using eight ensemble members and various synoptic background maps. 
Differences found in the predicted arrival time of each individual simulation were mostly due 
to CME initial speed and the time at which the CME was inserted at the WSA-ENLIL inner 
boundary, resulting in propagation through a different background solar wind. They used Na¬ 
tional Solar Observatory Global Oscillation Network Group (GONG) ( [Harvey et al, 199^ synoptic 
magnetograms and Air Force Data Assimilative Photospheric flux Transport (ADAPT) maps 
( jArge et al., 2010}|Henney et al., 2012 ). For their GME event they show that when using ADAPT 
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maps, the WSA-ENLIL model values were in better agreement with in-situ observations, and the ar¬ 
rival time predictions were improved due to the more accurate background solar wind representation. 
However, the overall spread in CME arrival times did not change significantly. 

This paper describes the WSA-ENLIL-bCone ensemble modeling system installed at the Com¬ 
munity Coordinated Modeling Center (CCMC) and results from the past 1.5 years of real-time 
execution at the CCMC/Space Weather Research Center. This is the first ensemble space weather 
prediction system for CME propagation of its kind employed in a real-time environment. The 
current version of the system evaluates the sensitivity of CME arrival time predictions from the 
WSA-ENLIL-bCone model to initial CME parameters. The CCMC, located at NASA Goddard 
Space Flight Center, is an interagency partnership to facilitate community research and accelerate 
implementation of progress in research into space weather operations. The SWRC is a CCMC 
sub-team which provides space weather services to NASA robotic mission operators and science 
campaigns, and prototypes new models, forecasting techniques and procedures. The CCMC also 
serves the CME Scoreboard websitc0 to the research community who may submit CME arrival time 
predictions in real-time for a variety of forecasting methods. The website facilitates model validation 
under real-time conditions and enables collaboration. For every CME event table on the site, the 
average of all submitted forecasts is automatically computed, thus itself providing a world-wide 
ensemble mean CME arrival time forecast from a variety of models/methods. 

In Section [2] a brief description of the WSA-ENLIL-bCone model is given. The triangulation 
algorithm for determining CME parameters for the ENLIL model is described in Section [3l The 
real-time ensemble modeling methodology is explained in Section |4] followed by an example of an 
ensemble simulation given in Section[5] Results and the evaluation of the first 1.5 years of simulations 
are described in Section [6l In Section [7] we discuss a parametric event case study of the sensitivity 
of the CME arrival time prediction to model free parameters for the CME and ambient solar wind. 
Finally, a summary and discussion are presented in Section [8l 


2. WSA-ENLIL+Cone Model Description 


The global 3D MHD WSA-ENLIL model provides a time-dependent description of the background 
solar wind plasma and magnetic field into which a CME can be inserted (Odstrcil, Smith, and Dryer, 1996H 


Odstrcil and Pizzo, 1999a||Odstrcil and Pizzo, 1999b||Odstrcil, 2003[[0dstrcil, Riley, and Zhao, 2004D .| 

This modeling system does not simulate CME initiation but uses kinematic properties of CMEs 
inferred from coronagraphs to launch a CME-like hydrodynamic structure into the solar wind 
and interplanetary magnetic field computed from the WSA coronal model (|Arge and Pizzo, 2000 


|Arge et ai, 2004 1. A common method to estimate the 3D CME kinematic and geometric param¬ 
eters is to assume that the geometrical CME properties are approximated by the Cone model 
( |Zhao, Plunkett, and Liu, 2002[ |Xie, Ofman, and Lawrence, 2004[ ) which assumes isotropic expan¬ 
sion, radial propagation, and constant CME cone angular width. Generally, a CME disturbance is 
inserted in the WSA-ENLIL model as slices of a homogeneous spherical plasma cloud with uniform 
velocity, density, and temperature as a time-dependent inner boundary condition at 21.5 solar radii 
{Rq) with an unchanged background magnetic field. While the simplest geometrical case is employed 
in this work, the WSA-ENLIL+Cone model can also support an elliptical geometry including tilt, 
an elongated spheroid or ellipsoid, and leading and trailing edge velocities. Measurements derived 


^ http: //ka uai.ccmc.gsfc.nasa.gov/CMEscoreboard 
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from coronagraphs (described in Section [3d1) determine the cloud velocity, location, and width. The 
CME cloud density (dcid) is a free parameter which by default is 4 times larger than typical mean 
values in the ambient fast wind providing a pressure of four times larger than that in the ambient 
fast wind. The cloud temperature is taken to be equal to the ambient fast wind temperature. 
Another ENLIL free CME parameter is the cavity ratio which allows the CME to be represented by 
a spherical shell of plasma and is based on coronagraph observations of CME cavities. The cavity 
ratio radcav is defined as the ratio of the radial CME cavity width to the CME width, with the 
default being no cavity (radcav=0). 

WSA-ENLIL+Cone runs performed for research and operations have shown that accurate de¬ 
scriptions of the heliosphere and transients are achieved only when the background solar wind is well- 
reproduced and if coronagraph observations from multiple views, for example from the SOlar and 
Heliospheric Observatory (SOHO) spacecraft near the Earth ( [Domingo, Fleck, and Poland, 1995| ) 
and the Solar TErrestrial RElations Observatory (STEREO) spacecraft ( Kaiser et al., 2008 1, are 
used to derive CME parameters ( [Lee et al., 2013|[Millward et al., 20T^ . WSA coronal maps provide 
the magnetic field and solar wind speed at the boundary between the coronal and heliospheric mod¬ 
els, usually at 21.5 Rq, and they are generated from synoptic magnetograms. Small latitudinal shifts 
in the magnetogram-derived coronal maps caused by inaccuracies in solar magnetic field observa¬ 
tions, particularly in the polar regions, can cause large longitudinal shifts in the solar wind structure, 
for example in characterizing high speed stream arrival times (e.g. IMacNeicel (I2009|l : IJian et ai\ 
([201ip b Other coronal models, such as MAS (MHD around a Sphere) ( Riley, Linker, and Mikic, 2001D [ 
or heliospheric tomography from interplanetary scintillation (IPS) (Jackson et al., 20111 can also 
provide the background solar wind and have been coupled with ENLIL heliospheric simulations. 

CCMC/SWRC has been carrying out routine WSA-ENLIL-|-Cone simulations for several years 
using solar magnetic synoptic maps and CME geometric and kinematic properties inferred from 
coronagraph observations ( jZheng et al., 2013 ). Each ENLIL run uses a WSA model synoptic map 
computed from the single GONG daily-updated synoptic magnetogram (see e.g. [Arge and Pizzo 


([20001) 1 closest to the time the simulation is executed. These low 4° resolution real-time simulations 
complete in ~20 minutes running on 2 nodes with 16 processors/node on a spherical grid size of 
256x30x90 (r, 9, </>) with 5-10 minute output cadence at locations of interest. The simulation range 
is 0.1 to 2 AU in radius r, -60° to 4-60° in latitude 6, and 0° to 360° in longitude tp. CME parameters 
are derived using real-time coronagraph observations from spacecraft and a geometric triangulation 
algorithm. The measurements are an approximation of the true 3D speed and width of the CME 
at 21.5i?Q (ENLIL inner boundary). However, often the coronagraph derived measurements are 
inferred from just a few data points, and some CMEs may be missed due to real-time data gaps. 
CME parameters derived in real-time and simulation graphical outputs are publicly available from 
the CCMC Space Weather Database Of Notifications, Knowledge, Information (DONKI) daXshas^ 


3. Ensemble CME Parameters 

3.1. STEREOCAT TRIANGULATION ALGORITHM FOR DETERMINING GME 
PARAMETERS 

CME parameters are determined using the Stereoscopic CME Analysis Tool (StereoCAT), developed 
by the CCMC for real-time CME analysis carried out by the CCMC/SWRC forecasting team. The 


^ http: //kauai.ccmc.gsfc.nasa.gov/DONKI 
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goal was to develop a tool that can be used quickly, yet reliably in a real-time environment with any 
possible combination of spacecraft available for analysis. It was also required that the tool was intu¬ 
itive and simple enough to be employed by a wide variety of users such as space weather forecasters, 
scientists, students, and citizen scientists. The basic methodology of the tool, i.e., tracking of CME 
kinematic properties from two different fields-of-view, is similar to that of the NOAA Space Weather 
Prediction Center CME Analysis Tool (CAT) bv IMillward et aD(|2013l) and the geometric localiza¬ 
tion developed bv IPizzo and Bieseckerl(|2004|l . However, StereoCAT does not attempt to capture the 
volumetric structure of CMEs but is based on tracking specific CME features. The algorithm is most 
similar to the CME geometric triangulation method of ILiu et'oTl (I2010|l . For a more detailed discus¬ 
sion of different CME analysis techniques in the context of cone model-based CME simulations, see 
dPulkkinen, Oates, and Taktakishvili, 20101 [Millward et ai, 2013D. Other stereoscopic methods for 


determining the kinematic properties of CMEs include those by Thernisien, Howard, and Vourlidail 

(|2006|) , [Lugaz et aT. d2010|) . and IDavies et aLI(l2013|) . 

StereoCAT is based on triangulation of transient CME features from two different coronagraph 
fields-of-view or planes-of-sky. We will call these planes-of-sky A and H, which may designate, for 
example, fields-of-view of the SECCHI/COR2 instruments onboard the STEREO A and STEREO 
B spacecraft ( Howard et ai, 2008[ ). The tool is used to manually identify the same CME features 
in two consecutive images which are then used to calculate the plane-of-sky velocities for A and H, 
and Vg, respectively. Note that these velocities are in local plane-of-sky coordinates indicated 
by ' and ". These data need to be brought into the same coordinate system (Heliospheric Earth 
Equatorial (HEEQ) coordinates in this case), which can be accomplished by rotations: 


VA = Ra ■ Va (1) 

v_B = Rb • (2) 

where operators Ra and Rs carry out transformations from planes-of-sky A and B coordinates 
into a common base such as HEEQ, respectively. 

We then define two projection matrices as 

Pa = 1 - eAe^ (3) 

Pb = 1 - eseg (4) 

where 1 is a 3 x 3 identity matrix. The unit vectors normal to the planes-of-sky of coronagraphs 
A and B are defined as ba and bb, where is the transpose of matrix ba. The matrices Pa and 
Pb project any vector to plane-of-skies A and B, respectively. Consequently, plane-of-sky speeds 
can be expressed as 


VA = Pa ■ V (5) 

vb = Pb ■ V (6) 

where v is the three dimensional vector pointing toward the propagation direction of the CME. 
While individual projection matrices are not invertible, we can combine Eqs. and (11]) to obtain 


(Pa + Pb) • v = va -k vb 


(7) 
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from which we can solve 


v=(Pa + Pb) ^-(va+vb) (8) 

Importantly, (Pa + Pb)-^ exists as long as planes-of-sky A and B are different, i.e. when ca and 
cb are not co-linear (parallel to each other). Therefore large triangulation errors occur when the 
spacecraft separation angle is very small or around 180°. 

A similar approach can be used to track the three-dimensional location r of a feature from 
plane-of-sky measurements ta and rg as 


r = (PA-fPB) ^-(rA-brs) (9) 

Often the time stamps of coronagraph imagery from spacecraft A and B do not match exactly. 
This is handled in StereoCAT by propagating the tracked feature in A with speed va to a new ta 
that matches the B time stamp. Consequently, matching time stamps are used for ta and in 
Eq. (ini). 

The angular size of a CME is estimated in StereoCAT simply by manually selecting the two 
outer edges of the CME. These two lines that connect through the center of the Sun are then used 
to compute the opening angle of the CME. It is noted that this process does not take into account 
projection of the outer CME edges to the spacecraft plane-of-sky, and is therefore a measurement 
of the projected CME width. While this is not an issue if the CME propagation direction is not 
too far away from the plane-of-sky of the spacecraft which is used to measure the opening angle, 
one needs to be very careful with events with propagation directions substantially away from the 
plane-of-sky, as in such cases the opening angle can be overestimated. This issue will be addressed 
in the future versions of StereoCAT. 

Other limitations of StereoCAT arise from the user’s ability to reliably identify the same struc¬ 
tures in images from both spacecraft due to ambiguities from the different viewing angles. It 
may at times be difficult or impossible to track the same structure since different sections of the 
CME contribute most strongly to images in different planes-of-sky ( [Howard and DeForest, 201 2D . 
Consequently, StereoCAT is not suitable to use with coronagraph data in which the CME appears 
as a halo, since the CME leading edge is not visible. 

3.2. PERFORMING CME MEASUREMENTS WITH STEREOCAT 


StereoCAT has three modes: “two-timepoint”, “ensemble”, and “frame series”, and is available 
online via a web interfaced ( jLaSota, 2013| ). Available coronagraphs include the LASCO C2 and 
C3 instruments on board the SOHO spacecraft (Brueckner et al., 19951, and the SECCHI/COR2 
instruments on the STEREO A and B spacecraft. All three modes are based on the same trian¬ 
gulation algorithm, described in section [3.II In the basic “two-timepoint” mode the user manually 
measures the CME leading edge height for two different times in each coronagraph image for two 
different coronagraph viewpoints. The plane-of-sky sky speed for each viewpoint is calculated, from 
which the triangulated speed and direction is computed using the algorithm described in section 
3.11 The user also manually measures the CME opening angle in each coronagraph view. Because 
this is a projected width measurement, both widths and their average are displayed for the user. 


^ http: //ccmc.gsfc.nasa.gov/analysis/stereo 
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Figure 1. Example screenshot of m=6 “two-timepoint” measurements performed for the 18 April 2014 CME using 
StereoCAT in “ensemble” mode. Two image pairs are shown from the SECCHI/COR2 instruments from STEREO 
B (a-b; top row) and STEREO A (c-d; bottom row), tor two different time steps, 2014-04-18 13:54 UT (a, c; left 
column) and 2014-04-18 14:24 UT (b, d; right column). The white circles indicate the 6 individual “two-timepoint” 
plane-of-sky leading edge height measurements (near the center of the CME front) and the width measurements are 
marked by the green circles (near the CME edges). The green lines in panel c illustrate the CME opening angle 
measurements for one of the coronagraph images. The plane-of-sky leading edge measurements (central white circles) 
are later combined together using the triangulation algorithm discussed in Sections 13.113.21 to generate 6^ = 36 
ensemble members. The distribution of the resulting CME parameters which are used as initial conditions for 36 
WSA-ENLIL-I-Cone simulations is shown in Figure |3] 


In “ensemble” mode the user manually repeats the same procedure as for the “two-timepoint” 
mode, by measuring the same feature for the same pair of coronagraphs at two different times. 
Between each “two-timepoint” measurement, the display is fully reset such that the user is forced 
to carefully remeasure the CME leading edge height and opening angle. This series of repeated 
measurements leads to a range of CME parameters which can be used to initialize an ensemble 
simulation. Eor every m “two-timepoint” measurements made, n = ensemble CME parameter 
members are automatically generated by combining different spacecraft measurement pairs. Eor 
example, for m = 2 “two-timepoint” measurements, there are n = 2^ = 4 ways to combine the first 
and second time step height measurements in viewpoints A and B to triangulate the CME. Since 
the two projected width measurements made for each measurement m are not triangulated, they 
are randomly assigned to each ensemble member. An example screenshot of m=6 “two-timepoint” 
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measurements performed in “ensemble” mode using StereoCAT is shown in Figure [H Two image 
pairs are shown from the SECCHI/COR2 instruments from STEREO B (top row) and STEREO 
A (bottom row), for two different times separated by 30 minutes in the left and right columns. The 
white circles indicate the 6 individual “two-timepoint” plane-of-sky leading edge height measure¬ 
ments (near the center of the CME front) and the width measurements are marked by the green 
circles (near the CME edges). The green lines in panel c of Figure [T] illustrate the CME opening angle 
measurements for one of the coronograph images. In this example the 6 individual “two-timepoint” 
measurements were combined by the algorithm to create 6^ = 36 ensemble members. 

After completing the measurements the user may inspect histograms of their CME parameters. 
The web interface allows the user to remove any ensemble members, and add any “custom” members. 
Generally, members are removed when they have nearly identical parameters, or members for which 
the triangulation appears unreliable. Custom members can be measurements from different image 
time pairs, from plane of sky estimates which incorporate the source location, or from any other 
CME measurement technique. The same procedure can be applied to create n individual ensemble 
measurements for x CMEs for a series of events which are then combined one-to-one to be simulated 
together such that there are n ensemble members containing x CMEs each. 

In “frame series” mode the user can measure a series of different frames (times) for each space¬ 
craft, which are then triangulated to create a CME height-time profile. The user selects a range 
of time, and steps through the images available from each instrument, measuring the CME in as 
many images as they choose. The software chooses time pairs of measurements for triangulation 
based on a user-specified maximum allowed time difference. From these measurements, plane-of-sky 
and triangulated height-time, velocity, latitude, and longitude profiles of the CME are generated. 
Triangulations made with different spacecraft pairs are shown as separate height-time profiles. 
Several methods are used to calculate the CME speed, acceleration, the time the CME passes 
21.5 Rs (ENLIL inner boundary), and the time it erupts from the Sun. These include least-squares 
linear and quadratic fits, averages over selected data points, and averages from only the first and last 
data points. Results for each method are reported separately, allowing the user to choose the most 
appropriate fitting technique depending on the acceleration profile of the CME. Plane-of-sky values 
are also reported, which can be used when coronagraph projection effects make this triangulation 
method unreliable. This can occur if the CME is very wide, appears as a halo, or is heavily projected 
in the coronagraph data. In these cases the user will not be able to identify the same CME leading 
edge feature in the data from two coronagraphs. The user can inspect the triangulated height values 
directly on the height-time plot to evaluate triangulation accuracy in these cases. 


4. Ensemble Modeling -with WSA-ENLIL-|-Cone 

The current implementation of this ensemble modeling method evaluates the sensitivity of WSA- 
ENLIL-|-Cone model simulations of CME propagation to initial CME parameters. As described in 
Section [3T] StereoCAT is used to create an ensemble of n CME parameters which are used as input 
to n WSA-ENLIL-I-Cone simulations. We have observed that n ~36 to 48 provides an adequate 
spread of input parameters, but this number can be increased if necessary. For n=48 a typical run 
takes 130 minutes to complete on 24 nodes with 4 processors/node on the initial development 
system. We estimate that the same run will take ~80 minutes on the CCMC production system 
that has 16 processors/node. 

The simulations provide n profiles of MHD quantities (density, velocity, temperature, and mag¬ 
netic field components) and a distribution of n predicted arrival times at locations of interest 
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within the computational domain. Currently, ensemble modeling is performed for spacecraft at 
the following locations: Mercury (MESSENGER), Venus (VEX), Earth (ACE, Wind, SOHO, and 
orbiting spacecraft). Mars (MSL, MAVEN, MEX), Spitzer Space Telescope, STEREO-A and B. The 
CME-associated disturbance/shock arrival time is then automatically computed in post-processing 
from any sharp increases in the modeled solar wind dynamic pressure at a given location. In this 
work, we focus on the ensemble results of the Earth-directed events. 

Eor Earth-directed CMEs, the CCMC/SWRC also computes n estimates of the geomagnetic Kp 
index using the WSA-ENLIL-|-Cone model plasma parameters at Earth. The geomagnetic three- 
hour planetary K index, Kp, is a measure of general planetary wide geomagnetic disturbances 
at mid-latitudes based on ground-based magnetic observations ( [Bartels, Heck, and Johnston, 193^ 
[Rostoker, 1972[ [Menvielle and Berthelier, 1991| ). The Kp index is created from standardized K 
indices from individual stations, which measure the magnitude of horizontal geomagnetic field 
disturbances (not including daily variations). Kp is a quasi-logarithmic index ranging from 0 to 
9. Real-time estimated planetary Kp indices are available from NOAA using real-time data from 
a limited number of geomagnetic observatories, and the final definitive Kp is from the Helmholtz 
Center Potsdam GEZ German Research Centre for Geosciences. 

The predicted Kp estimate is made by using the INewell et al. \ (|2007|) coupling function arising 
from their correlation of 20 candidate coupling functions with geomagnetic indices. The function 
which represents the rate of magnetic flux d^Mp/dt opening at the magnetopause and correlated 
best with 9 out of 10 indices is given as 

d^Mp/dt = ■!;buik'‘'^^-BT'^^sin®/^(^), (10) 

where fbuik is the bulk solar wind speed, the interplanetary magnetic field (IME) clock angle dc 
is given by tan~^ [By/B^), and the perpendicular component of the magnetic field is given by 
Bt = [By + (in GSM coordinates). An exponential fit to the correlation of this coupling 

function with the A'p index yields the following relation used for the estimate 

_ g g _ g2.17676—5.2001(c!<I>Mp/iii) (11) 

lEmmons et al.\ (|2013p showed for their sample of 15 events, that Kp predictions using Eq. |TT] 
computed directly from in-situ solar wind observations had a mean absolute error of 0.5. Because 
ENLIL modeled CMEs do not contain an internal magnetic field and the magnetic field amplihcation 
is caused mostly by plasma compression, only the magnetic field magnitude is used and three 
magnetic field clock angles scenarios of 90° (westward), 135°(south-westward), 180° (southward) 
are assumed. This provides a simplistic estimate of three possible maximum values which the Kp 
index might reach following arrival of the predicted CME shock/sheath. Eor the forecast, the Kp 
estimates are rounded to the nearest whole number. 

Another commonly used activity index is the Dst (disturbance storm time) index, which is a 
measure of magnetosphere storm activity primarily from the strength of the ring current. The 
index is obtained from the measurement of the perturbations in the horizontal component of 
the Earth’s magnetic field from ground-based observatories that are sufficiently distant from the 
auroral and equatorial electrojets, are located at approximately ±20° geomagnetic latitude, and 
are evenly distributed in longitude ( jSugiura, 1964| ). Although the ring current makes the largest 
contribution to the Dst, all magnetospheric current systems contribute, such as the Chapman- 
Ferraro magnetopause current which is strengthened during sudden storm commencement (SSC) 
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Figure 2. Coronagraph observations of the 18 April 2014 CME with an onset time at 13:09 UT as viewed from (b) 
SOHO LASCO/C2 and C3, (a) STEREO SECCHI/COR2 B and (c) A near the time of 14:50 UT. The fields-of-view 
of LASCO/C2, C3, and SECCHI/COR2 are 2.2-6i?Q, 2.8-32 /?q (shown here cropped to 17i?0), and 2.5-15R0 
respectively. Images from helioviewer (http://www.helioviewer.org) ||Muller et al., 2009 |l. 


and increases the Earth’s surface field and gives a sudden positive jump in Dst. Currently, ENLIL 
model results are not used to predict the Dst, however in principle this can be computed in a similar 
manner to the Kp index by using the INewell et al.\ (j2007|l Dst relation. 


5. Example Ensemble: 18 April 2014 CME 

In this section we describe the real-time ensemble modeling of an Earth-directed partial halo CME 
that was first observed at 13:09 UT on 18 April 2014 by by SECCHI/COR2A. Figure [5] shows 
this CME as viewed from SOHO LASCO/C2 and C3, STEREO SECCHI/COR2 A, and B near 
14:50 UT. This CME was associated with an M7.3 class solar flare from Active Region (AR) 12036 
located at S18°W29° with peak at 13:03 UT. The eruption and a coronal wave were visible south 
of the active region in SDO/AIA 193A and a nearby filament eruption was visible in AIA 304A. 
Subsequently starting at 13:35 UT, an increase in solar energetic particle proton flux above 0.1 
pfu/MeV (1 pfu = 1 particle cm“^ sr“^ s“^) was observed by the GOES-13 EPEAD instrument in 
Earth orbit. 

Figure [T] shows StereoCAT measurements for the 18 April 2014 CME. As discussed above, the 
central white circles indicate the individual leading edge measurements and the green outer circles 
near the CME edges are the projected width measurements. The six leading edge measurements 
are combined together using the triangulation algorithm discussed in Sections 13.1113.21 to generate 
6^ = 36 ensemble members. The distribution of the resulting CME parameters which are used 
as initial conditions for n=36 WSA-ENLIL-|-Cone simulations is shown in Figure [3] in the (a) 
equatorial plane (latitude=0°) and (b) meridional plane (longitude=0°). The plots show the CME 
velocity vectors in spherical HEEQ coordinates with the grids showing the degrees longitude (a) 
and latitude (b), and the radial coordinate showing the speed in km/s. The Sun-Earth line is along 
0° longitude and latitude. The arrow directions on the grid indicate the CME central longitude 
and latitude respectively, with CME half width indicated by the color of the vector. The arrow 
lengths correspond to the CME speed. CME propagation directions are clustered between -30 to 
-40° latitude, and around 10° west of the Sun-Earth line in longitude, while CME speeds range from 
^1300 to 1600 km/s. Median CME parameters are: speed of 1394 km/s, direction of 9° longitude, 
-35° latitude, and a half-width of 46°. 
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Figure 3. Distribution of the 18 April 2014 CME input parameters shown in the (a) equatorial plane (latitude=0°) 
and (b) meridional plane (longitude=0°). The plots show the CME speed vectors in spherical HEEQ coordinates 
with the grids showing the degrees longitude (a) and latitude (b), and the radial coordinate showing the speed in 
km/s. The Sun-Earth line is along 0° longitude and latitude. The arrow directions on the grid indicate the CME 
central longitude and latitude respectively, with CME half width indicated by the color of the vector. The arrow 
lengths correspond to the CME speed. CME propagation directions are clustered between -30 to -40° latitude, and 
around 10° west of the Sun-Earth line in longitude, while CME speeds range from ^^1300 to 1600 km/s. Median 
CME parameters are: speed of 1394 km/s, direction of 9° longitude, -35° latitude, and a half-width of 46°. 
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Figure 4. Global view of the 18 April 2014 CME on 20 April at 06:00 UT: WSA-ENLIL+Cone scaled velocity 
contour plot for the (a) constant Earth latitude plane, (b) meridional plane of Earth, and (c) 1 AU sphere in 
cylindrical projection, for the ensemble member with median CME input parameters (speed of 1394 km/s, direction 
of 9° longitude, -35° latitude, and a half-width of 46°). Panel (d) shows the measured (red) and simulated (blue) 
radial velocity profiles at Earth, with the simulated CME duration shown in yellow. 
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WSA-ENLIL+Cone Ensemble Profiles ond ACE Observotions ot Eorth 



Figure 5. 18 April 2014 CME ensemble: Model calculated density, velocity, magnetic field, and temperature profiles 
at Earth for all 36 ensemble members plotted as color traces along with the observed in-situ LI observations from 
ACE plotted in black (red for Bz). The model traces are color coded by CME input speed such that slow to faster 
input speeds are colored from light green to dark blue. The observations show clear signatures of the arrival of 
an ICME, including a leading shock (abrupt increase in all the solar wind parameters at around 10:20 UT) with 
enhanced post-shock temperatures, enhanced magnetic field with rotations in direction, and declining solar wind 
speed. The spread in the color traces show that most of the predictions are earlier than the observed arrival, with a 
mean predicted arrival at Earth of 20 April 2014 at 05:07 UT and a range from 20 April 2014 at 01:08 UT to 11:16 
UT. 
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Ensemble modeling of CMEs using the WSA-ENLIL+Cone model 


Model results for the 36-meinber ensemble WSA-ENLIL+Cone run for this CME are shown in 
EiguresIMSl Eor the ensemble member with median CME input parameters, Eigure|4] shows a scaled 
velocity contour plot for the (a) constant Earth latitude plane, (b) meridional plane of Earth, and 
(c) 1 AU sphere in cylindrical projection on 20 April at 06:00 UT. Panel (d) shows the measured 
(red) and simulated (blue) radial velocity profiles at Earth, with the simulated CME duration 
shown in yellow. This simulation figure shows the northeastern portion of the CME impacting 
Earth. Eigure [5] shows the modeled magnetic field, velocity, density, and temperature profiles at 
Earth plotted as color traces for all 36 ensemble members, along with the observed in-situ LI 
observations from ACE, plotted in black. The model traces are color coded by CME input speed 
such that slow to faster input speeds are colored from light green to dark blue. The arrival of the 
CME-associated shock was observed by Wind and ACE on 20 April 2014 at around 10:20 UT, 
and energetic storm particles were observed by ACE. The provisional SYM-H index (^1 minute 
Dst) shows a sudden storm commencement of +25 nT at 11:01 UT. The observations in Eigure[5] 
show clear signatures of the arrival of an Interplanetary Coronal Mass Ejection (ICME), including a 
leading shock (abrupt increase in all the solar wind parameters at around 10:20 UT) with enhanced 
post-shock temperatures, enhanced magnetic field with rotations in direction, and declining solar 
wind speed. This CME was predicted to arrive at Earth and also at Mars for all of the 36 runs. The 
mean predicted arrival at Earth was on 20 April 2014 at 05:07 UT with arrival times from individual 
runs ranging from 20 April 2014 at 01:08 to 11:16 UT. A histogram showing the distribution of 
arrival times at Earth is shown in Eigure [S] with individual arrivals marked by the blue arrows. 
This figure shows a normal distribution with 50% of the predicted arrivals within one hour of 
the mean. The prediction error for the mean predicted CME arrival time was -5.2 hours and the 
observed arrival time was just within the ensemble predicted spread. The spread in ensemble member 
predictions can also be seen in Figure [5] compared to the observations, showing that most of the 
predictions are earlier than the observed arrival with a few after. From the CME input parameters 
plotted in Figure[3]the ensemble members with arrival times closest to the observed time had CME 
input speeds in the range of 1200-1400 km/s, latitudes near -40° and half widths around 35°-40°. 
This suggests that the early arrival time predictions for this event could be due to overestimations 
of the CME input speed and half width. 

The NOAA real-time observed Kp index (and the Potsdam final Kp) reached 5 during the 
synoptic period 12:00-15:00 UT on 20 April associated with the CME shock arrival. The Dst reached 
a minimum of -24 nT at 15:00 UT on 21 April and thus based on Dst, this CME only resulted 
in very weak geomagnetic activity. As discussed in Section SI Eq. [TT] can be used to forecast the 
maximum Kp index from maximum ENLIL predicted quantities at CME shock/sheath arrival at 
Earth (colored traces shown in Figure [5]). Figure [7] shows the predicted probability distribution 
of Kp for three clock angle scenarios 6c = 90° (green), 135° (purple), 180° (orange). The figure 
also shows the overall Kp forecast probability distribution calculated for all three angles combined 
90°-180°, assuming each scenario is equally likely, in black. The standard deviation of the overall 
Kp forecast probability distribution is 1.1, with 84% of the forecasts falling between Kp = 5 to 
7. The most likely forecast is for Kp=7 at 41%, followed by Kp=5 at 27% and Kp=6 at 16% 
likelihood of occurrence. Using the most likely forecast of Kp=7, the Kp prediction error for this 
event is AKpei-r = Kppj-edicted ~ Kpodserved = 2 (overprediction). The overprediction oi Kp may be 
related to the overestimation of the CME input speed. In Sections I6.HI6.2I and [7] we discuss various 
factors which can contribute to early arrival time predictions and Kp overpredictions. 


SOLA: MaysEnsembleAccepted.tex; 13 May 2015; 1:03; p. 13 





Mays et al. 


Ensemble of ENLIL CME orrivols at Eorth 



Averoge Arrivol: 201A-04-20 05:07UT 


Figure 6. 18 April 2014 CME: Histogram distribution of arrival time predictions at Earth (bin size of 1 hour) with 
individual arrivals marked by the blue arrows. This figure shows a normal distribution with 50% of the predicted 
arrivals are within one hour of the mean. The prediction error for the mean predicted CME arrival time is -5.2 hours 
and the observed arrival time was within the ensemble predicted spread. 


6. Real-time Ensemble Modeling: First Results 

For 35 Earth-directed CME events from January 2013 through June 2014, real-time ensemble 
modeling was carried out by the CCMC/SWRC team following the methods described in Sections 
[41151 In Table [T] we list a summary of the ensemble simulation results for these 35 CME events. The 
first and second columns give the CME onset date and time based on the first appearance in C2 or 
COR2. Generally, if two CMEs occur within a day of each other they will both be included in the 
same simulation as separate CMEs which may or may not merge during their propagation. A few 
of the ensemble simulations listed in the table contain two CMEs as part of a single run. In these 
cases, CMEs that were simulated together with the CME listed on the previous row are indicated by 
*. The third column lists (for 2013) the 2nd order plane-of-the-sky (POS) speed at 20 i?© reported 
in the SOHO LASCO CDAW CME catalog ( [Yashiro et al., 2004[ [Gopalswamy et ai, 2009[ ). If 
measurements were not made to 20 Rq, the 2nd order POS speed at the time of last observation 
is used. The next four columns provide the median ensemble CME input parameters of v, latitude, 
longitude (HEEQ), and half-width {w/2) measured using StereoCAT. In columns 8, 9 and 10 we 
list the mean predicted arrival time of all ntot ensemble members, followed by the spread in arrival 
times in hours relative to the mean. The next column (11) shows Upredicted hits; the number of 
ensemble members out of ntoti the total number of ensemble members, that predict that the CME 
will arrive at Earth. This ratio p = Upredicted hits/»^tot gives a forecast probability and conveys 
the forecast uncertainty about the likelihood that the CME will arrive. Columns 12, 13, and 14, 
list the actual arrival time of the CME-associated shock or disturbance observed in-situ at the 


^http: //cdaw.gsfc.nasa.gov/CME_list 
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Ensemble modeling of CMEs using the WSA-ENLIL+Cone model 


Probabilistic Kp Forecast Distribution (18 April 2014) 
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Figure 7. Distribution of Kp probability forecast using ENLIL predicted solar wind quantities at Earth for three 
clock angle scenarios Oq = 90° (green), 135° (purple), 180° (orange), and all three angles combined 90°-180° (black) 
(assuming equal likelihood). The standard deviation of the overall Kp forecast probability distribution is 1.1, with 
84% of the forecasts falling between Kp = 5 to 7. The most likely forecast is for Kp=7 at 41%, followed by Kp=5 
at 27% and Kp=() at 16% likelihood of occurrence. The NOAA real-time observed Kp index (and the Potsdam final 
Kp) reached 5 during the synoptic period 12:00-15:00 UT on 20 April associated with the CME disturbance arrival. 



4 .5 6 7 8 9 
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Wind spacecraft, followed by the total in-situ observed CME transit time relative to the CME 
start time. In the last column the prediction error Aterr is calculated for predictions indicating 
hits. The prediction error is defined as Aterr = ^predicted — ^observed, which is negative when ENLIL 
predictions are earlier than the observed CME arrival time, and late predictions are positive. When 
possible, ICME and magnetic cloud catalogs were used to help assess whether the CME did arrive 
at Earth. These included the IRichardson and Canel (I20I0|1 ICME catalog, and the Wind ICME 
catalofUwith circular flux rope model fitting (based on |Hidalgo et aL| ()2000|l b Shocks identified by 
the SOHO CELIAS/MTOE/PM “shockspotter” program were also used in arrival time assessment. 
Determining the measured in-situ arrival time of the CME-associated shock or disturbance can be 
subjective and therefore can be a source of error in the prediction error calculation. Taking this 
into consideration, in-situ signatures which could not be unambiguously identified as the arrival of 
the CME-related disturbance are indicated by f and these five ensembles are not included in the 
following forecast verification. This reduces the sample size from 35 to 30 ensembles. 

In the following subsections we discuss ensemble CME arrival and Kp forecast verification 
inspired by methods used in ensemble weather forecasting and applied here for the hrst time. 


® http: //www. srl.caltech.edu/ACE/ASC/DATA/level3/icmeta ble2.htm 
® http: //wind. nasa.gov/index_WLI CM EJist.htm 
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6.1. CME ARRIVAL FORECAST VERIFICATION 


To begin with a simple forecast evaluation of CME arrival time, the ensemble mean can be taken 
as a single forecast. Using the prediction error Aterr = ^predicted — ^observed (last column of Table [U 
the mean absolute error (MAE) is 12.3 hours, the Root Mean Square Error (RMSE) is 13.9 hours, 
and the mean error (ME) is -5.8 hours (early) for all 17 ensembles containing hits. Considering the 
sample size in this study, these errors are comparable to CME arrival time prediction errors (a RMSE 
of ^ 10 hours) reported by others ( Millward et al., 2013[ Romano et al., 2013[|Vrsnak et ai, 2014[ 
[Mays et al., 2014[ ). Similarly, [Colaninno, Vourlidas, and Wu (|2013l) used a variety of methods to 
evaluate CME arrival time predictions (non real-time) based on imaging data analysis only, and 
found an error ±6 hours for 78% of their sample, and ±13 hours for their full sample of 9 CMEs. 
The CME arrival time prediction error is inevitably related to the CME propagation speed, thus 
it is useful to consider the input speed and in-situ observed transit time relative to the prediction 
error. For this sample, the average in-situ observed transit time was 66 hours. In Figure [HK the 
CME arrival time prediction error is plotted against the CME input speed, and in Figure [SJd the 
prediction error as a percentage of the CME transit time is plotted against the CME input speed. 
The error bars are computed using the predicted ensemble range as listed in column 10 of Table 
[I] The dashed horizontal line indicates the mean arrival time prediction error (a) and mean of 
the prediction error/transit time percentage (b). These figures show a nearly consistent negative 
prediction error for fast CMEs above ~1000 km/s such that these fast CMEs are generally predicted 
to arrive earlier than they are observed. This could be a sign of the modeled CME having too much 
momentum as defined by a combination of the input speed and half width (which is related to 
the modeled CME mass). The overestimation of the modeled CME velocity compared to in-situ 
observed values is also due to the modeled CME having a lower magnetic pressure than is observed 
in typical magnetic clouds. 
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Table 1.: Summary of the ensemble simulation results for 35 CME events (January 2013 - June 
2014). Columns 1-2: CME onset date and time. Column 3: SOHO LASCO CME Catalog plane-of- 
sky (POS) speed for 2013. Columns 4-7: median ensemble CME input parameters of v, latitude, 
longitude (HEEQ), and half-width (w/2). Columns 8-10: mean predicted arrival time of all ritot 
ensemble members, and the spread in arrival times in hours relative to the mean. Column 11: 
?^predicted hits, the number of ensemble members predicting that the CME will arrive at Earth out of 
ntot, the total number of ensemble members. Columns 12-14: actual arrival time observed in-situ, 
and the observed CME transit time relative to the CME start time. Column 15: prediction error 
Aterr = ^predicted ~ ^observed for hits, Or CR and FA for correct rejections and false alarms. 


CME Onset 

SOHO 

Median CME parameters 

Mean Predicted Arrival 

Date Time 

^POS 

V 

Lat 

Lon 

to/2 

Date 

Time ! 

Spread 

(yyyy-mm-dd) 

(UT) 

(km/s) 

(km/s) 

(°) 

n 

n 

(UT) 


(h) 

2013-01-13 

07:24 

229t 

342 

1 

10 

28 

2013-01-17 

06:30 

+C>.b 

-6.2 

2013-01-16 

19:00 

616 

750 

-26 

52 

42 

2013-01-19 

21:33 

+ 11.1 
-16.3 

2013-02-06 

00:24 

1851 

1460 

30 

-26 

30 

2013-02-08 

05:37 

+8.7 

-6.8 

2013-04-11 

07:24 

819 

1000 

0 

-15 

55 

2013-04-13 

06:14 

+6.1 

-5.4 

2013-06-21 

03:12 

1903 

1997 

-15 

-48 

60 

2013-06-22 

13:02 

+5.7 

-3.6 

2013-06-30 

02:24 

410t 

386 

9 

4 

34 

2013-07-02 

20:56 

+ 1.1 
-0.6 

2013-07-16 

04:00 

639 

795 

-19 

9 

19 

2013-07-18 

20:29 

+6.9 

-6.2 

2013-08-02 

13:24 

443 

596 

-16 

28 

13 



2013-08-07 

18:24 

473 

570 

-25 

11 

44 

2013-08-11 

05:03 

+6.8 

-5.9 

2013-08-08 

23:54 

41lt 

454 

-17 

14 

18 



2013-08-30 

02:48 

884 

861 

21 

-48 

59 

2013-09-01 

08:34 

+4.6 

-5.4 

2013-09-19 

03:36 

4491 

362 

-15 

-43 

27 



2013-09-29 

20:40 

1164 

1000 

26 

30 

66 

2013-10-02 

04:11 

+9.1 

-10.8 

2013-10-06 

14:39 

710t 

747 

1 

2 

16 

2013-10-09 

22:10 

+ 10.1 
-9.8 

2013-10-22 

04:24 

697 

764 

51 

-10 

49 

2013-10-25 

08:19 

+ 10.2 
-10.9 

2013-12-04 

2013-12-05* 

23:12 

00:00 

585 

623+ 

697 

651 

41 

25 

-9 

63 

46 

31 

2013-12-07 

13:45 

+4.1 

-4.1 

2013-12-12 

2013-12-12* 

03:36 

06:24 

943 

695 

1067 

694 

-32 

-52 

51 

8 

50 

50 

2013-12-14 

18:11 

+ 16.1 
-13.9 

2013-12-29 

2013-12-29* 

00:12 

05:48 

296+ 

477 

682 

495 

32 

-33 

8 

-58 

47 

43 

2014-01-01 

02:39 

+ 13.2 
-9.3 


Pi = 
^hits/ 

In-situ Arrival Transit 

Date Time Time 

(UT) (h) 

Aterr 

(h) 

20/48 

2013-01-16 23:25+ .... 


18/48 

2013-01-19 16:48 69.8 

4.8 

19/48 

2013-02-08 03:15+ 


36/36 

2013-04-13 22:13 62.8 

-16.0 

47/48 

2013-06-23 03:51 48.7 

-14.8 

4/36 


CR 

28/48 

2013-07-18 12:55+ 


0/24 


CR 

48/48 


FA 

0/48 


CR 

46/48 

2013-09-02 01:56 71.1 

-17.4 

0/24 


CR 

36/36 

2013-10-02 01:15 52.6 

2.9 

22/24 

2013-10-08 19:40 53.0 

26.5 

45/47 


FA 

37/48 


FA 

36/48 

2013-12-15 16:30+ 


48/48 

48/48 


FA 


* CMEs was simulated together with the CME listed on the previous row as part of a single ensemble, 
t 2nd-order plane-of-sky speed at last possible measured height. 

I In-situ signature could not be unambiguously identified as arrival of CME-related disturbance and is not included in forecast verification. 
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Table 1.: Continued from previous page 


CME Onset 

SOHO 

Median CME parameters 

Mean Predicted Arrival 

Pi = 

In-situ Arrival 

Transit 


Date Time 

^^POS 

V 

Lat 

Lon 

io/2 

Date 

Time Spread 

'^hits/ 

Date 

Time 

Time 

^terr 

(yyyy-mm-dd) 

(UT) 

(km/s) 

(km/s) 

(°) 

(°) 

n 

(UT) 


(h) 

^7tot 

(UT) 


(h) 

(h) 

2014-01-07 

18:24 

1714 

2399 

-28 

38 

64 

2014-01-09 

00:17 

-ty.2 

-6.9 

48/48 

2014-01-09 

19:39 

49.3 

-19.4 

2014-01-30 

16:24 

940 

843 

-50 

-28 

45 

2014-02-02 

10:10 

+ 11.9 
-12.3 

13/24 

2014-02-02 

23:20 

78.9 

-13.2 

2014-01-31 

15:39 

370 

718 

12 

-29 

40 

2014-02-03 

17:20 

+9.1 

-6.0 

+22.8 

-18.03 

12/12 

23/24 




FA 

2014-02-04 

2014-02-04* 

01:09 

16:24 

501 

323 

778 

560 

-35 

-34 

21 

23 

49 

42 

2014-02-06 

20:39 

2014-02-07 

16:28* 


2014-02-12 

05:39 

494 

740 

8 

5 

59 

2014-02-14 

23:47 

+ 13.2 
-8.7 

48/48 

2014-02-15 

12:46 

79.1 

-13.0 

2014-02-18 

01:25 

712 

882 

-24 

-43 

52 

2014-02-20 

16:29 

+ 17.7 
-28.7 

29/36 

2014-02-20 

02:42 

49.3 

13.8 

2014-02-19 

16:00 

510 

883 

-32 

-10 

29 

2014-02-22 

12:20 

+ 13.2 
-12.6 

32/36 

2014-02-23 

06:09 

86.2 

-17.8 

2014-02-25 

01:09 

2069 

1394 

-18 

-80 

78 

2014-02-26 

22:15 

+20.7 

-11.3 

40/48 

2014-02-27 

15:50 

62.7 

-17.6 

2014-03-23 

03:36 

834 

715 

-5 

-60 

55 

2014-03-26 

00:58 

+ 11.8 
-12.6 

38/48 

2014-03-25 

19:10 

63.4 

5.80 

2014-03-23 

06:12 

536 

503 

37 

-45 

34 



0/12 

2/36 

14/16 




OR 

2014-03-29 

18:12 

505 

707 

36 

41 

43 

2014-04-01 

21:30 

+1.0 

-1.0 

+6.3 

-10.1 




CR 

2014-04-02 

13:36 

1527 

19 

-55 

51 

2014-04-04 

19:01 

2014-04-05 

10:00 

68.4 

-15.0 

2014-04-18 

13:09 


1394 

-35 

9 

46 

2014-04-20 

05:07 

+6.1 

-4.0 

36/36 

2014-04-20 

10:20 

45.2 

-5.2 

2014-06-04 

15:48 


580 

-40 

-28 

50 

2014-06-07 

20:56 

+6.8 

-8.1 

22/36 

2014-06-07 

16:12 

72.4 

4.7 

2014-06-10 

13:09 


980 

-9 

-89 

64 

2014-06-12 

20:28 

+3.1 

2/36 

12/12 




CR 

2014-06-19 

2014-06-30 

17:12 

07:24 


569 

751 

3 

-12 

-20 

-63 

44 

29 

2014-06-22 

2014-07-02 

16:12 

20:32 

—.1 
+7.6 
-5.4 
+7.4 
-12.3 

2014-06-22 

18:28 

73.3 

-2.3 

CR 

0/36 






* CMEs was simulated together with the CME listed on the previous row as part of a single ensemble. 

I In-situ signature could not be unambiguously identified as arrival of CME-related disturbance and is not included in forecast verification. 
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Ensemble modeling of CMEs using the WSA-ENLIL+Cone model 


(a) CME arrival time prediction error compared to CME input speed 


(b) Prediction error relative to CME transit time and input speed 
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Figure 8. (a) CME arrival time prediction error plotted against the CME input speed, (b) Prediction error as a 
percentage of the CME transit time, plotted against the CME input speed. The error bars are computed using the 
predicted ensemble range as listed in Tabled 


Table 2. Forecast performance contingency table for 30 
ensembles. 



CME arrival forecast 

Observation 

Will occur 

Will not occur 

Occurs 

Hit (17) 

Miss (0) 

Does not occur 

False alarm (5) 

Correct rejection (8) 


Ensemble modeling produces a probabilistic forecast pi of the likelihood of CME arrival for 
each ensemble but we begin with a more simple forecast evaluation by binning the probability 
Pi into a categorical yes/no forecast. Categorical forecasts only have two probabilities, zero and 
one. Therefore we start by binning the probability forecast p into two categories: “yes” the CME 
will arrive, and “no” the CME will not arrive. In the signal detection theory model of weather 
forecasting, event forecasting performance can be evaluated in terms of a 2x2 contingency table, 
as shown in Table [5] ( [Harvey et al., 199^ [Weigel et al, 2006) [Jolliffe and Stephenson, 201 1[ ). Eor 
CME arrival prediction, the “event” is taken as the “CME arrival”. Hits are then defined as CME 
arrivals which were both predicted and observed to occur. Misses are defined as CME arrivals 
which were not predicted, but were observed to occur. False alarms (FA) are defined as CME 
arrivals that were predicted to occur, but were observed not to occur. And correct rejections (CR) 
are CME arrivals that were not predicted, and were observed not to occur. To bin each ensemble’s 
probabilistic forecast, correct rejections were identified when the criterion of the forecast probability 
Pi = ?^predicted hits/'^-totai members < 15% was met; i.e., that less than 15% of the total predictions in 
the ensemble indicated CME arrival. Similarly, the inverse criterion is used to identify hits. Table 
[2] shows the contingency table definitions and values for this 30 event sample: 17 hits, 8 correct 
rejections, 5 false alarms, and 0 misses (see Tabled) for specific CR and FA events). For this sample 
zero misses indicates that there were no ensemble simulations which did not predict CME arrivals 
which were observed to occur. There were 8 out of 30 correct rejections and 5 false alarms for events 
that were not observed in-situ, giving a correct rejection and false-alarm rate of 62% (8/ 13) and 
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38% (5/ 13) respectively. The correct alarm ratio, defined as the number of hits over the number 
of hits and false alarms, is 77% and the false alarm ratio is 23%. 

Let us now consider a more nuanced technique to evaluate the probabilistic forecast with¬ 
out partitioning it into a categorical forecast with only two probabilities as described above. A 
method defining the magnitude of probability forecast errors is the Brier Score (BS) ( |Brier, 1950[ 
[Murphy, 1973| [Wilks, 1995| ), defined as 


N 




( 12 ) 


where N is the number of events, pi is the forecast probability of occurrence for event i, and is 1 
if the event was observed to occur and 0 if it did not occur. For CME arrival prediction, the “event” 
here is taken as the “CME arrival” and pi is listed in column 11 of Tabled] for each ensemble. This 
score is a probability mean square error which weights larger errors more than small ones and ranges 
from 0 to 1, with 0 being a perfect forecast. The BS computed from all iV= 30 ensemble CME 
arrival probabilities (Tabled! column 11) is 0.15, which indicates that in this sample, on average, 
the probability p of the CME arriving is fairly accurate. However, such verification scores reduce 
the problem to a single measure which can only consider one dimension, whereas there are many 
dimensions to the system. Eor example, consider the aspect of forecast reliability. Reliable forecasts 
are those where the observed frequencies of events are in agreement with the forecast probabilities. 
To evaluate the reliability of probabilistic ensemble forecasts, a set of probabilistic forecasts pi 
must be evaluated using observations that demonstrate that those events either occurred or did 
not occur. Multiple forecasts must be evaluated because a single probabilistic forecast cannot be 
simply assessed as “right” or “wrong” e.g. if a forecast suggests a 30% chance of CME arrival, and 
the CME does arrive, the forecast is not clearly either “right” or “wrong”. Therefore, to provide 
forecast verification for a p=30% chance of CME arrival one would need to compile the statistics 
of observed CME arrivals for a set of forecasts that predicted a 30% chance of arrival. In this way 
a reliability diagram can be constructed to determine how well the predicted probabilities of an 


event correspond to their observed frequencies ( [Wilks, 1995[ [ Jolliffe and Stephenson, 2011[ ). Figure 
[n^ shows the reliability diagram of the likelihood of CME arrival forecast for the 30 event sample, 
with the reliability for this sample shown as the black line with points and the perfect reliability 
diagonal as a dotted line. The line of perfect reliability is diagonal because, for example, when a 60% 
probability forecast is made, it is considered perfectly reliable if the event is observed to occur 60% 
of the time over multiple ensemble forecasts. The number of events used in each calculation is shown 
next to each point, and the sample size is smaller than needed for a robust diagram. Nevertheless, 
the diagram shows that overall ensemble modeling is underforecasting in the forecast bins between 
20-80%, and slightly overforecasting in the 1-20% and 80-100% forecast bins. Overforecasting is 
when the forecast chance of CME arrival (forecast probability) is higher than is actually observed; 
i.e., the CME is observed to arrive less often than is predicted. Similarly, underforecasting is when 
the chance of CME arrival is lower than is actually observed; i.e., the CME is observed to arrive 
more often than is predicted. 

Another aspect of forecast reliability is to assess how well the ensemble spread of the fore¬ 
cast represents the true variability of the observations. For 8 out of 17 of the ensemble runs 
containing hits the observed CME arrival was within the spread of ensemble arrival time pre¬ 
dictions. This indicates that roughly half of the observations fall outside of the extremes of the 
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Figure 9. CME arrival time forecast verification: (a) Reliability diagram of the forecast probability of CME arrival 
for the 30 ensemble sample, with the ensemble results shown as the black line with points and the diagonal perfect 
reliability as a dotted black line. The number of ensembles used in each calculation is shown next to each point. 
The diagram indicates underforecasting in the forecast bins between 20-80%, and slight overforecasting in the 1-20% 
and 80-100% forecast bins. Overforecasting is when the forecast probability of CME arrival is higher than observed; 
i.e. the CME is observed to arrive less often than is predicted. Similarly, underforecasting is when the CME arrival 
forecast probability is lower than observed; i.e. the CME is observed to arrive more often than is predicted, (b) Rank 
histogram for the 17 ensembles containing hits indicates undervariability of initial conditions. 


predicted ensemble spread. However, one aspect of a reliable forecast is that the set of ensem¬ 
ble member forecast values for a given event and observations should be considered as random 
samples from the same probability distribution. This reliability then implies that if an n member 
ensemble and the observation are sorted from earliest to latest arrival times, the observation is 
equally likely to occur in each of the n -I- 1 possible “ranks”. Therefore a histogram of the rank 
of the observation, “rank histogram”, tallied over many events should show be uniform (flat) 
( [Anderson, 1996[ [Hamill and Colucci, 1997| [Talagrand, Vautard, and Strauss, 1997D . While more 
samples would be desirable, it is still instructive to examine the rank histogram for the CME 
arrival time predictions from the 17 ensembles containing hits in this sample, shown in Figure [Hb- 
Since each ensemble run in our sample does not have the same number of members, the rank has 
been normalized to 10 (9 member ensemble). To construct this rank histogram the CME arrival time 
predictions of each ensemble are sorted from earliest to latest and the rank of where the observed 
arrival falls among the predicted times is noted. For example, an ensemble with a rank of 8 has the 
meaning that 7 arrival time predictions fall before the observed arrival, a rank of 10 would mean 
that all 9 predictions occur before the observation, and a rank of 1 means that the observation 
occurs before all of the predictions. The non-uniform U-shape of this histogram partly illustrates 
that roughly half of the observed arrivals are outside the spread of predictions (ranks 1 and 10), 
with a tendency for an overall early spread of predictions (rank=10) compared to observations 
(also quantified by mean arrival time error of -7.0 hours). U-shaped rank histograms can indicate 
lack of variability in the ensemble, but can also be a sign of a combination of conditional biases in 
the model ( [Hamill, 2001| ). However, when evaluating the WSA-ENLIL-bCone model in this sample 
of ensembles, and >70 regular runs containing hits performed by SWRC ([Romano et al., 2013[), 
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an overall negative bias (early predictions) was found, with less bias for CME input speeds below 
~1000 km/s. Therefore, it is unlikely that a combination of positive and negative model biases 
within the ensembles contributed to the U-shaped rank histogram for our sample. Most likely, the 
U-shape suggests undervariability, indicating that these ensembles to not sample a wide enough 
spread in CME input parameters. 

6.2. Kp FORECAST VERIFICATION 

For each event for which a hit is predicted in Table [T] ensemble modeling provides a probabilistic 
Kp forecast (see Section |4|) for three magnetic field clock angles scenarios of 90° (westward), 135° 
(south-westward), 180° (southward). An overall probabilistic Kp forecast can then be obtained 
by making the simple assumption that each clock angle is equally likely to occur. Table [3] lists 
the overall probabilistic Kp forecast p(Kp = h) for each Kp bin b (e.g. the distribution shown in 
Figure[7]in black) for these 17 events. The observed Kp, sudden storm commencement (SSC) and 
minimum Dst indices are also shown. The mean predicted Kp is listed in column 12, along with the 
overall predicted Kp spread (using plus or minus notation). Underlined Kp probabilities indicate 
that the NOAA real-time observation falls within this bin, and the final definitive Kp values are 
listed in column 13. The Dst values are from the real-time (quicklook) Dst index provided by the 
World Data Center for Geomagnetism in Kyoto, Japan. In order to estimate the reliability of the 
probabilistic Kp forecast, the Brier Score is calculated for each Kp bin and listed on the last line 
of the table. 

To evaluate forecast performance, a single categorical predicted Kp forecast can be derived from 
the probabilistic Kp forecast p{Kp = b) distribution. For example, the single categorical forecast 
^'^Ppredicted Can be taken as the mean predicted Kp, or the most probable Kp value. This allows 
a Kp prediction error to be computed as AKp^,.,. = Kpp;.g(ji(.tg(j — Kpobsgj-^gj for each ensemble, 
where positive values of AKp^^^ indicate an over prediction of the Kp index and negative values 
indicate that Kp has been under predicted. If the categorical Kppj.g(jigtg(j is taken as the Kp bin 
b which has the highest likelihood in the probabilistic Kp forecast p{Kp = b) for each ensemble, 
the prediction errors are calculated to give a mean absolute error (MAE) of 1.9, Root Mean Square 
Error (RMSE) of 2.5, and mean error (ME) of -1-1.4. However, if the categorical Kppj.g(jigtg(j is taken 
as the mean predicted ATp in each ensemble (last column of Table [3]) these errors are reduced to 
MAE=1.5, RMSE=2.0, and ME=-|-0.6. Consequently, utilizing the ensemble mean Kp yields a 
more accurate forecast in this sample, however both forecast choices show an overall tendency for 
the overprediction of Kp. Given that the modeled CMEs do not have an internal magnetic field 
structure, the INewell et al.\ (|2007l) Kp coupling function using ENLIL results as input performs 
surprisingly well. For comparison, using ACE solar wind data as input to the coupling function for 
this sample gives Kp prediction errors of MAE=0.67, RMSE=0.77, and ME=-|-0.22. 

In FigureflTlthe Kp prediction error (from the ensemble mean Kp) is compared to the CME input 
speed; the error bars show the ensemble Kp prediction spread. This figure shows that Kp is usually 
overpredicted when CME input speeds are above ^1000 km/s. This bias is also apparent in the Kp 
predictions made from a sample of >70 regular WSA-ENLIL-|-Cone runs reported bv IRomano et aO 
(|2013|) . The Kp overprediction is most likely due to an overestimation of the CME dynamic pressure 
at Earth by the WSA-ENLILj-Cone model, due to the CME having a lower magnetic pressure than 
is observed in typical magnetic clouds. Also since the CME dynamic pressure is linearly related 
to the density and the square of the velocity, this quantity will be in particular more sensitive to 
higher CME input speeds, and produce higher in-situ speeds than those measured. Another factor 
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Table 3. Summary of Kp prediction results for 17 ensemble runs containing hits. Columns 1-2: CME start date 
and time. Columns 3-11: overall probabilistic Kp forecast for each Kp bin assuming equal likelihood of three 
clock angle scenarios. Underlined Kp probabilities indicate that the NOAA real-time Kp observation falls in 
this bin and the observed definitive Kp is in column 13. The mean predicted Kp is listed in column 12, along 
with the overall predicted Kp spread (using plus or minus notation). The Brier Score BS is calculated for each 
Kp bin and listed on the last line of the table. The Dst sudden storm commencement and minimum values are 
listed in the last two columns. 


CME Onset 


Binned Probabilistic Kp Forecast (%) 


Mean Kp 

Obs. 

Dst (nT) 

Date Time (UT) 

1 

2 

3 

4 

5 

6 

7 

8 

9 

S>c spread 

Kp 

ssc 

min. 

2013-01-16 

19:00 

0 

13 

26 

28 

6 

11 

9 

4 

4 


4- 

+6 

-34 

2013-04-11 

07:24 

0 

0 

0 

0 

0 

33 

5 

62 

0 


3+ 

+21 

-7 

2013-06-21 

03:12 

0 

0 

0 

0 

4 

16 

23 

43 

15 

7 +2 
‘ -2 

5+ 


-49 

2013-08-30 

02:48 

0 

0 

6 

31 

28 

33 

2 

0 

0 


3+ 


-31 

2013-09-29 

20:40 

0 

0 

6 

26 

24 

39 

5 

0 

0 

St? 

8- 

-1-30 

-67 

2013-10-06 

14:39 

0 

67 

33 

0 

0 

0 

0 

0 

0 

2 

^ -0 

6- 

+21 

-65 

2014-01-07 

18:24 

0 

0 

0 

6 

8 

19 

25 

26 

16 

7 +2 

' -3 

3- 

+2 

-14 

2014-01-30 

16:24 

0 

0 

15 

13 

33 

13 

18 

8 

0 


2+ 

+15 

-7 

2014-02-12 

05:39 

0 

0 

12 

25 

40 

24 

0 

0 

0 


5o 

-1-52 

-16 

2014-02-18 

01:25 

0 

1 

10 

21 

29 

26 

10 

2 

0 


6o 


-86 

2014-02-19 

16:00 

0 

2 

30 

34 

28 

5 

0 

0 

0 

4^2 

4+ 

+4 

-56 

2014-02-25 

01:09 

0 

0 

1 

11 

16 

21 

22 

21 

9 


5+ 


-99 

2014-03-23 

03:48 

0 

0 

16 

28 

28 

24 

4 

0 

0 

4t? 

4- 

-1-20 

-18 

2014-04-02 

13:36 

0 

0 

21 

19 

40 

12 

7 

0 

0 

4t? 

4o 

+16 

-16 

2014-04-18 

13:09 

0 

0 

0 

3 

27 

17 

40 

14 

0 


5o 

+25 

-24 

2014-06-04 

15:48 

0 

0 

18 

29 

36 

17 

0 

0 

0 

4t? 

6+ 

+31 

-38 

2014-06-19 

17:12 

0 

0 

8 

31 

33 

25 

3 

0 

0 

4t? 

3o 

+14 

-9 

Brier Score 

0.00 

0.09 

0.27 

0.19 

0.17 

0.17 

0.02 

0.04 

0.00 




in the higher CME dynamic pressure can arise from the approximation of the CME as a cloud with 
homogeneous density in the model. 

Other factors contributing to Kp overprediction may include the magnetic field direction - 
two out of the three field configurations assumed produce persistent southward fields (135° and 
180°), so there is a bias towards geoeffective field configurations. Examining the distribution of 
north-south magnetic fields associated with the ICMEs of IRichardson and Canel (I2010p and the 
associated sheaths, in only 2% of cases are southward fields completely absent, so the bias towards 
geoeffective field configurations is consistent with observations. However, both small and large 
maximum southward fields are observed relatively infrequently (e.g., maximum southward fields 
are <4 nT in 17% of events, and >15 nT in 16%), suggesting that the weighting of 90° and 
180° clock angles should be reduced. In particular, reducing the 180° clock angle weight would be 
expected to reduce the Kp overprediction. 

The last line of Table [3] lists Brier Score calculated for each Kp bin. Here, the BS is a measure 
of the magnitude of error in the Kp probability forecast (how likely a given Kp bin will occur) 
in each bin. The BS values indicate that in this sample, the Kp probability forecast is reliable 
for the Kp=5 and 6 bins {BS=0.17 for both), and less so for the Kp=3 and 4 bins {BS=0.27 
and 0.19). Although the scores also indicate that the forecast is most reliable for the smallest and 
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Figure 10. Kp forecast verification: (a) Histogram of the observed Kp values (black) and the forecast Kp 
probability distribution (hashed) for this sample (see Tablc[3ll. (b) Rank histogram of Kp predictions for all ensembles. 
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Figure 11. Kp prediction error (computed from the ensemble mean Kp) compared to the CME input speed shows 
an overprediction of the Kp value for CME input speeds above ~1000 km/s. Error bars indicate the ensemble Kp 
prediction spread listed in [3] 


largest Kp bins, most of the observations in this sample did not fall in these extreme bins, hence 
a larger sample is needed to verify forecast reliability for these bins. Figure ITOh shows the overall 
observed Kp distribution and the forecast Kp probability distribution for the events in Table [3] 
used to calculate the BS. 

To further evaluate Kp probability forecast reliability we compare the observed Kp to the 
spread in ensemble predictions. For most (12 out of 17) of the ensembles, the observed Kp was 
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Figure 12. Coronagraph observations of the 11 April 2013 07:24 UT CME as viewed from (b) SOHO LASCO, (a) 
STEREO SECCHI/COR2 B and (c) A near the time of 09:55 UT. 


within the overall predicted Kp spread (column 11). The observed Kp was also within the predicted 
mean Kp±l for 11 out of 17 of the ensembles. A rank histogram was also constructed for the Kp 
predictions for all ensembles and is shown in Figure [TOji, again normalized to an ensemble size of 
9. To construct this rank histogram the Kp predictions are sorted from smallest to largest and the 
rank of where the observed Kp value falls among the predicted Kp values is noted. For example, 
an ensemble with a rank of 6 has the meaning that 5 Kp predictions are less than the observed 
Kp value, a rank of 10 would mean that all 9 of the Kp predictions are less than the observed Kp 
(underprediction), and a rank of 1 means that the observed Kp value is less than all of the Kp 
predictions (overprediction). The histogram has an overall flat shape, with more occurrence at rank 
1 (the observed Kp was less than the predicted range) and less occurrence in the higher ranks which 
shows the bias for Kp overprediction (mean error=+0.6). Note, that the rank histogram does not 
indicate how “good” forecasts are but only measures whether the observed probability distribution 
is well represented by the ensemble. Therefore, a uniform, flat rank histogram is a necessary but 
not sufficient condition for determining the reliability of ensembles ([Hamill, 2001[). 


7. ENLIL Parameter Sensitivity: 11 April 2013 Event Case Study 

In the current configuration, other than the measured CME speed, direction, and size, the real-time 
WSA-ENLIL-I-Cone ensemble simulations use the default values for the model CME free parameters. 
In this section, we present a case study which examines the effect of changing these model free 
parameters on the ensemble modeling. The CME starting on 11 April 2013 at 07:24 UT was chosen 
for this study due to the large early arrival time prediction error obtained for all members of the 
model ensemble. |Taktakishvili, MacNeice, and Odstrcil| (j2010|l studied the dependence of arrival 
time predictions on the uncertainty in CME input parameters (speed, width, density ratio) for 
three Earth directed CME events of varying speeds. A similar procedure was adopted for this case 
study, and by employing the ensemble modeling technique, the parameter space can be sampled 
more systematically. 

The original set of simulations performed in real-time were chosen as the “base ensemble” 
(ensemble I). Subsequently, ten ensemble runs (ensembles II-XI), each containing 36 members for 
360 total simulations, were performed to assess the sensitivity of the CME arrival time prediction 
to changes in the model free parameters and ambient solar wind model, while keeping the CME 
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Figure 13. Distribution of the 11 April 2013 CME input parameters shown as speed vectors (all CME half widths are 
55°), in the same format as Figure[3] Median CME parameters are: speed of 1000 km/s, direction of -15° longitude, 
0° latitude, and a half-width of 55°. This figure shows that custom ensemble members were chosen with speeds of 
850, 900, 1000, 1100, and 1200 km/s, between ±-10° latitude, -10° to -25° longitude with a half width of 55°. 


speed and direction input parameters fixed. The ENLIL model free parameters considered in this 
study include the CME half width, CME density ratio, CME cavity ratio, and ambient solar wind 
reduction factor. The CME density ratio (dcid) is a free parameter which is by default is a set factor 
of 4 times larger than typical mean values in the ambient fast wind providing a pressure of four 
times larger than that of the ambient fast wind. The cavity ratio radcav is defined as the ratio of 
the radial CME cavity width to the CME width, with the default being no cavity radcav=0. The 
ambient speed reduction factor vred reduces the solar wind speed provided by the WSA coronal 
map in order to account for expansion of the solar wind from the WSA boundary to 1 AU since 
WSA is calibrated against 1 AU in-situ observations. 

Figure [12] shows the CME starting on 11 April 2013 at 07:24 UT as viewed from SOHO LASCO 
C3, STEREO A and B SECCH1/COR2 near the time of 09:55 UT. On this date the STEREO 
B spacecraft was located at -142° and STEREO A was at 133° in HEEQ coordinates. This CME 
was associated with an M6.5 class flare from AR 11719 located at N07E13 with peak intensity at 
07:16 UT. The eruption, coronal dimming and wave were visible mostly southeast of the active 
region in SDO/AIA 193A. Additionally, an increase in solar energetic particle proton flux was 
observed starting at around 07:40 UT by the SOHO COSTEP (reaching 1 pfu/MeV, in the 16- 
40 MeV energy range), ACE EPAM (100 pfu/MeV, 1.22-4.94 MeV), and GOES-13 EPEAD (5 
pfu/MeV, 15-40 MeV energy range) instruments starting at 08:00 UT, and by the IMPACT HET 
instruments on STEREO B (5 pfu/MeV, 24-41 MeV energy range) and A (0.001 pfu/MeV, 24-41 
MeV energy range). This solar energetic particle event and its longitudinal extent is studied in 
detail bv ICohen et a/.l(l2014|) and lLario et aZ.ip014l) . 

Due to the lack of availability of real-time concurrent coronagraph images, triangulation of CME 
parameters with the StereoCAT “ensemble” mode method was not possible for this CME. Therefore, 
the ensemble was composed of “custom” members. The CME parameters for each member were 
derived from plane-of-sky CME speed measurements combined with the source location at the 
Sun. The distribution CME input parameters for 36 ensemble members are shown in Figure 1131 
Median CME parameters are: speed of 1000 km/s, direction of -15° longitude, 0° latitude, and a 
half-width of 55°. This figure shows that custom ensemble members were chosen with speeds of 
850, 900, 1000, 1100, and 1200 km/s, between ±10° latitude, -10° to -25° longitude with a half 


SOLA: MaysEnsembleAccepted.tex; 13 May 2015; 1:03; p. 26 












Ensemble modeling of CMEs using the WSA-ENLIL+Cone model 


2013-04-13106:00 

(a) Constant latitude plane 
W90 


(b) Meridional plane 
N90 Lon = 0° 


2013-04-11T00 +2.25 days 

(c) Radial plane (lon-lat plot) 


R= 1.00 AU 
















' 





\ 

j 


r 





-1 













7 



> 














W90 W180 



Vr (km/s) f 


measured simulated 


200 550 900 1250 1600 - 

ENLIL-2.7 lowres-a3b11 WSA_V2.2 GONG-2135 


Figure 14. Global view of 11 April 2013 CME on 13 April at 06:00 UT: WSA-ENLIL+Cone scaled velocity contour 
plot in the same format as[4|for the ensemble member with median CME input parameters (speed of 1394 km/s, 
direction of 9° longitude, -35° latitude, and a half-width of 46°). 


width of 55°. Subsequent re-analysis of the CME height-time evolution gives average plane-of-sky 
speeds of ~800 km/s and ~700 km/s for SECCHI COR2B and LASCO C3 respectively, yielding 
a triangulated speed of 850±200 km/s, -5°±5° latitude, -15°±10° longitude, 50°±5° half width, 
which is represented within the ensemble members derived in real-time. 

The WSA-ENLIL-|-Cone model scaled velocity contour plot is shown in Eigure [TT] on 13 April 
at 06:00 UT for the ensemble member with median CME input parameters. This simulation figure 
shows a nearly direct CME impact at Earth, slightly eastward. Figure [12] shows the base ensemble 
WSA-ENLIL-|-Cone modeled quantities for all 36 ensemble members (color traces) at Earth along 
with in-situ ACE (black) and Wind (grey) observations (when there are ACE datagaps). The model 
traces are color coded by CME input speed such that slow to faster input speeds are colored from 
light green to dark blue. All 36 of the ensemble members predicted that the CME would arrive 
(100%) and the mean predicted arrival at Earth was 13 April 06:14 UT (range from 13 April 
00:47 to 12:20 UT). The histogram of the distribution of arrival times is shown in Figure ITOl The 
clustering of predicted arrival times in this histogram (and also in Figure [T21) reflects the limited 
number of discrete CME input speeds represented in the ensemble (see Figure ESI), with faster 
CMEs arriving first. The CME-associated shock was observed to arrive at ACE and Wind on 13 
April at 22:13 UT, giving an average prediction error of-16 hours. Clear ICME signatures including 
enhanced low variability magnetic field, declining solar wind speed, and low proton temperatures, 
start at around 16:45 UT on 14 April through about 18:30 UT on 15 April. The overall spread 
in arrival time predictions of all of the members in the base ensemble (including the clustering by 
CME input speed) can also be seen in Figure ESI as the color traces increase ahead of the observed 
arrival. The traces also show that the velocity, density, and temperature are overpredicted, while 
the maximum magnetic field strength is similar to that actually observed. The passage of this CME 
did not produce a geomagnetic storm due to an almost persistently northward magnetic field, shown 
in red in the top panel of Figure ESI The NOAA real-time observed Kp index reached 3 during 
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Figure 15. 11 April 2013 CME base ensemble: Model calculated magnetic field, velocity, density, and temperature 
profiles at Earth for each ensemble member along with the observed in-situ LI observations from ACE in black 
(red for Bz). Wind density observations are plotted in grey due to missing ACE values. The model traces are color 
coded by CME input speed such that slow to faster input speeds are colored from light green to dark blue. The 
CME-associated shock was observed to arrive at ACE and Wind on 13 April at around 22:13 UT with clear ICME 
signatures starting around 16:45 UT on 14 April through about 18:30 UT on 15 April. All of the arrival times 
indicated by the model results are earlier than the observed shock arrival, and are clustered by CME input speed. 
The mean predicted arrival time at Earth is 13 April 06:14 UT, with a range from 13 April 00:47 to 12:20 UT. 
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Ensemble of ENLIL CME arrivals at Earth 



Figure 16. 11 April 2013 base ensemble: Histogram of distribution of arrival time predictions at Earth (one hour 
bin size). The actual arrival was observed on 13 April at around 22:13 UT by Wind. The clustering of predicted 
arrival times reflects the limited number of different CME input speeds represented in the ensemble (see Figure fT^ . 
with faster CMEs arriving first. 


the synoptic period 21-24:00 UT on 13 April, while the Potsdam final Kp was 3-I-. The Dst index 
shows a sudden storm commencement of -1-21 nT at 23:00 UT on 13 April, and reached a minimum 
of only -7 nT at 11:00 UT on 15 April. 

In Figure flTl the arrival time prediction error (Aferr = ^predicted — ^observed) for the members in 
the base ensemble is plotted against the CME input speed for different CME input propagation 
directions (gray scale coded) and a fixed half width of 55° (full angular width of 110°). On 11 
April 2013 Earth was located at -5.9° latitude and 0° longitude in HEEQ coordinates, thus the 
input propagation direction of -10° latitude and -10° longitude (black, and dark blue in subsequent 
figures) represent the members with the most direct impact. This Figure shows that the arrival 
time prediction error ranges from -9.9 hours to -21.4 hours and increases with initial CME speed. 

Considering that a source of the prediction error could be due to uncertainty in the CME width, 
an identical ensemble (11) simulation was performed with the same input conditions but decreasing 
the half width by 10° (full angular width decreased from 110° to 90°). Figure [15] shows the difference 
from the original predicted arrival times from the base ensemble against the CME input speed for 
different propagation directions (as shown in Figure flTl) when the full angular width decreased from 
110° to 90°. In this figure (and those subsequent), the new CME arrival time prediction error is 
shown in hours relative to the original “base ensemble” prediction error. Compared to the original 
arrival time estimates, the overall prediction error decreases by 0.2 to 1.8 hours with increasing initial 
CME speed. Since all of the predictions in the base ensemble were too early (negative prediction 
error), a decrease in prediction error means that the new predictions are shifted to later times, closer 
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14 April 2013 Base Ensemble: Prediction Error vs. Input Speed for Different Input Propagation Directions 



Figure 17. 11 April 2013 base ensemble: CME arrival time prediction error (Aterr = ^predicted ~ ^observed) for tbe 
ensemble members plotted against the CME input speed for different CME input propagation directions (gray scale 
coded). 


to the observed arrival time. Nevertheless, the improvement is small compared to the prediction 
error. 

Next, the dependence of prediction on the input CME density ratio dcid was considered. Two 
ensembles were performed (III and IV) for which all parameters of the base ensemble were held 
fixed but the CME density ratio was adjusted from four (default), to two, and three. The results are 
shown in Figure [19] which shows the difference from the predicted arrival time for the base ensemble 
as a function of CME speed for the two different density ratios. The prediction error decreases by 
3.3 to 4.3 hrs for a CME density ratio of two and by 1.3 to 1.7 hrs for a density ratio of three, as 
a function of increasing initial CME speed. Hence, reducing the density ratio from four to two or 
three improves the arrival time prediction by around 3.5 or 1.5 hours, respectively. 

Another ENLIL CME parameter is the cavity ratio which allows the CME to be represented 
by spherical shell of plasma, based on coronagraph observations of CME cavities. The cavity ratio 
radcav is defined as the ratio of the radial CME cavity width to the CME width, with the default 
being no cavity radcav =0. Figure [201 shows the results of five ensembles (V-IX) with the CME cavity 
ratio adjusted to 0.1, 0.3, 0.5. 0.6, and 0.7, i.e., the CME is modeled as a progressively thinner shell 
as the ratio increases, using the base ensemble for all other parameters fixed. Specifically, the 
differences from the arrival times obtained for the base ensemble are plotted as a function of CME 
speed and direction (indicated by the symbol/line type) for each of these ensembles (indicated by 
the line color). For a cavity ratio of 0.1 the prediction remains largely unchanged compared to the 
base ensemble (with 0.15 hours). For the other cavity ratios, increasing differences in Figure [20l as 
the cavity ratio increases correspond to the predicted arrival moving to later times, reducing the 
prediction error. Furthermore, for each cavity ratio there is a spread (1 to 3 hours) in prediction time 
difference (compared to the base ensemble) for the different CME input directions with the more 
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14 April 2013 Ensemble II: Decrease Full Width by 20° 



Figure 18. 11 April 2013 ensemble II: Difference from the base ensemble (as shown in Figure fTTl l in hours when 
the CME input half width is decreased by 10° against the CME input speed for different propagation directions. 


14 April 2013 Ensembles III & IV: Decrease CME Density Ratio (default dcld=4) 



Figure 19. 11 April 2013 ensembles III and IV: Difference from the base ensemble (as shown in Figure fTTl in hours 
when the CME density ratio dcid is decreased to dcld=3 and dcld=2 (default dcld=4) against the CME input speed 
for different propagation directions 
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14 April 2013 Ensembles V-IX: Increase CME Cavity Ratio (default radcav=0) 



Figure 20. 11 April 2013 ensembles V-IX: Difference from the base ensemble (as shown in Figure fTTl) in hours when 
the CME cavity ratio is increased (radcav=0.1, 0.3, 0.5. 0.6, and 0.7) against the CME input speed for different 
propagation directions. Different CME input directions are indicated by the symbol/line type and each ensemble is 
indicated by a different line color. The cavity ratio radcav is defined as the ratio of the radial CME cavity width to 
the CME width, and the default is no cavity radcav =0. 


Earth directed inputs showing the largest difference from the base ensemble. Overall the prediction 
error decreases by between 0-1.6 hrs, 0.9-3.9 hrs, 1.8-4.7 hrs, 2.4-5.6 hrs for cavity ratios of 0.3, 0.5, 
0.5, 0.6 respectively. 

Considering now the influence of the ENLIL ambient solar wind solution, the in-situ data-model 
comparison (Figure [TS]) for the base ensemble indicates that the modeled background solar wind 
speed is '^125 km/s higher than the observed in-situ values, whereas the default value of ambient 
solar wind reduction factor vred is 25 km/s . To examine the role of the speed reduction factor in 
the prediction error, vred was increased to 50 km/s and 75 km/s for two ensembles (X and XI). 
This factor reduces the speed provided by the WSA coronal map in order to account for expansion 
of the solar wind from the WSA boundary to lAU since WSA is calibrated against 1 AU in-situ 
observations. Figure ITU shows the prediction time difference from the base ensemble for these two 
ensembles which show differences of 1-1.3 hrs and 2.2-2.8 hrs for vred of 50 km/s and 75 km/s 
respectively. Since the differences are positive, this indicates that the predicted arrival times are 
moved later, reducing the error relative to the observed arrival time. As might be expected, the 
modeled CME propagates more slowly when the ambient solar wind is slower. Figure [22] illustrates 
how the modeled background solar wind speed better matches the observed speed prior to CME 
arrival when the ambient speed reduction factor vred is increased to 75 km/s in ensemble XL 

Overall, this parametric case study shows that after the CME input speed, the cavity ratio and 
density ratio assumed in ENLIL have the greatest effects on the predicted CME arrival time, each 
changing this time by about 3 hours on average. Their effect is also more noticeable with higher 
CME input speeds. The CME input speed, cavity and density ratios define the CME momentum 
which defines the CME deceleration. In the addition to the using the default values, new ensemble 
runs could be performed with changes to the CME cavity ratio and density ratio as informed 
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14 April 2013 Ensembles X & XI; Increase Ambient Speed Reduction (default 25 km/s) 



Figure 21. 11 April 2013 ensembles X and XI: Difference from the base ensemble (as shown in Figure fTTll in hours 
when the ENLIL ambient speed reduction factor vred is increased to 50 km/s and 75 km/s (from the default value 
of 25 km/s) against the CME input speed for different propagation directions. 



Figure 22. 11 April 2013 ensemble XI: The predicted velocities (color traces) better match the observed in-situ 
values at ACE (black) when vred is increased to 75 km/s compared to the default vred =25 km/s results shown in 
Figure fTSi The second peaks in the predicted solar wind speed are artifacts of the ENLIL modeled CME as a spherical 
cloud. 
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by coronagraph measurements of the CME. Here, we have only examined the effect of changing 
the ad hoc ambient speed reduction factor in ENLIL, but we could also produce an ensemble of 
ambient solar wind WSA-ENLIL simulations using different ambient speed reduction factors which 
can be compared to in-situ measurements to determine “best” factor to use in subsequent CME 
simulations. An ensemble forecast reflecting uncertainties in the background solar wind could also 
be produced by using a variety of magnetograms (from different observatories or processed using 
different techniques) as input to the WSA or WSA-ADAPT models. 


8. Summary and Discussion 

This study evaluates the first ensemble CME prediction system of its kind employed in a real-time 
environment, providing unique space weather information for NASA users. The ensemble prediction 
approach provides a probabilistic forecast which includes an estimation of arrival time uncertainty 
from the spread in predictions and a forecast confidence in the likelihood of CME arrival. The 
current implementation explores the sensitivity of CME arrival time predictions from the WSA- 
ENLIL-|-Cone model to initial CME parameters. First results give a mean absolute arrival time error 
of 12.3 hours, RMSE of 13.9 hours, and mean error of -5.8 hours (early bias), based on a sample 
of 30 CME events for which ensemble simulations were performed. The arrival time is generally 
based on the arrival of the CME-generated shock at the Earth. The ensemble mean absolute error 
and RMSE are both comparable with other CME arrival time prediction errors reported in the 
literature. 

When considering the overall performance of CME arrival prediction, it was found that the 
correct rejection rate is 62%, the false-alarm rate is 38%, correct alarm ratio is 77%, and false 
alarm ratio is 23%. Each ensemble CME arrival time forecast includes a forecast probability p = 
J^predicted hits/^T-tot, which conveys a forecast uncertainty about the likelihood that the CME will 
arrive, which can be compared with observations to determine forecast reliability. The Brier Score 
[BS) of 0.15 for all 30 ensemble CME arrival probabilities indicates that in this sample, on average, 
the predicted probability of the CME arriving is fairly accurate. (A BS of 0 on a range of 0 to 1 is a 
perfect forecast.) However, the reliability diagram (Figure |9^) shows that the ensemble simulations 
are underforecasting the likelihood that the CME will arrive in the forecast bins between 20-80%, 
and slightly overforecasting in the 1-20% and 80-100% forecast bins. Overforecasting is when the 
forecast chance of CME arrival is higher than is actually observed; i.e., the CME is observed to 
arrive less often than is predicted. More ensemble simulations are needed for a more robust forecast 
verification of these probabilistic CME arrival time forecasts. 

For 8 out of 17 of the ensemble runs containing hits, the observed CME arrival was within the 
spread of ensemble arrival time predictions. The initial distribution of CME input parameters was 
shown to be an important influence on the accuracy of CME arrival time predictions. Particularly, 
the median and spread of the input distribution should accurately represent the range of CME 
parameters derived from observations. This is evidenced by the rank histogram (Figure [Sjo) which 
illustrates that roughly half of the observed arrivals are outside the spread of predictions, and also 
suggests undervariability in initial conditions; i.e., these ensembles do not sample a wide enough 
spread in CME input parameters. 

Each set of ensemble simulations also provides a probabilistic Kp forecast p(Kp = h) for each Kp 
bin b which can be compared with observations to determine forecast reliability. The Brier Score 
{BS) for the probabilistic Kp forecast bins show reliability for the Kp=5 and 6 bins {BS=0.17 for 
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both), and less so for the Kp=3 and 4 bins {BS=0.27 and 0.19). If choosing a single categorical 
Kp forecast value, the mean predicted Kp was found to have smaller prediction errors compared 
to using the Kp bin with the highest likelihood from the probabilistic Kp forecast. The observed 
Kp was within ±1 of the predicted mean Kp for 11 out of 17 of the ensembles. The Kp prediction 
errors computed from the mean predicted Kp show a mean absolute error of 1.4, RMSE of 1.8, 
and mean error +0.4. There is a known overall tendency for the overprediction of Kp, generally 
found for CME input speeds above 800-1000 km/s. Again, more ensemble simulations are needed 
to provide better forecast verification and to calibrate the Kp forecast. 

This paper focuses on the forecast verification of the ensemble modeling aspect of CME arrival 
and Kp predictions. More events, as well as comparison of results using different CME propagation 
models, are needed for more comprehensive forecast verification. These aspects are being investi¬ 
gated in a separate verification study which evaluates >400 single WSA-ENLIL+Cone simulations 
(of which there are >70 simulations containing CME arrivals) performed at the CCMC/SWRC. 

The parameter sensitivity studied discussed in Section [7] suggests future directions for this en¬ 
semble system. In the addition to the using the default model values, new ensemble runs could be 
performed with changes to the CME cavity ratio and density ratio as informed by coronagraph 
measurements of the CME. As discussed in Section [2] an accurate representation of the background 
solar wind is necessary for simulating transients, and prediction errors arising from background 
characterization and other model limitations should be considered. An ensemble forecast reflecting 
uncertainties in the background solar wind could be produced by using a variety of magnetograms 
(from different observatories or processed using different techniques) as input to the WSA or 
WSA-ADAPT models. From these results one can produce an ensemble of ambient solar wind 
WSA-ENLIL model outputs which can be compared to in-situ measurements to determine “best” 
coronal maps/model instance. These sub-selected WSA or WSA-ADAPT maps could then be used 
for a series of ensemble WSA-ENLIL+Cone CME simulations. Such an improved ensemble forecast 
would produce predictions which also reflect uncertainties in the WSA-ENLIL modeled background 
solar wind in addition to the uncertainties in CME input parameters (as considered in this work). 

Another improvement could be the use of real-time interplanetary scintillation (IPS) observations 
by the Ooty Radio Telescope ( Manoharan, 2006| ). These data can provide crucial information about 
the CME propagation and interaction with the surrounding solar wind which could be used to 
provide updated information on CME parameters as the CME moves out from the Sun. This 
information could then be used to refine model predictions of the propagation of the CME. The 
STEREO Heliospheric Imagers also provide CME propagation information out to lAU. However, it 
is not always possible to extract this information from real-time data, and the imagers do not always 
have an optimal viewing angle for Earth-directed CMEs. Comparisons of CME propagation from 
WSA-ENLIL with near real-time observations of the CME location inferred from IPS, the STEREO 
heliospheric imagers, or some other source, can be used to select ensemble members with the best 
agreement using quantitative and visual inspection employing advanced visualization techniques 
such as “3D volumetric rendering” ( |Bock et ai, 2014| ). 

Finally, the forecasting of CME arrival would benefit from the use of other propagation models, in 
addition to WSA-ENLIL, each with its own set of independently assessed input parameters, leading 
to a community-wide ensemble prediction capability. A first step to such a capability is provided 
by the CME Scoreboard, described in Section [I] where anyone is invited to post their estimate of 
the arrival time of a recently observed CME in real-time. 
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